Topology Of Gene Regulatory Networks Research Articles

Correlated stabilizing selection shapes the topology of gene regulatory networks.

The evolution of gene expression is constrained by the topology of gene regulatory networks, as co-expressed genes are likely to have their expressions affected together by mutations. Conversely, co-expression can also be an advantage when genes are under joint selection. Here, we assessed theoretically whether correlated selection (selection for a combination of traits) was able to affect the pattern of correlated gene expressions and the underlying gene regulatory networks. We ran individual-based simulations, applying a stabilizing correlated fitness function to three genetic architectures: a quantitative genetics (multilinear) model featuring epistasis and pleiotropy, a quantitative genetics model where each genes has an independent mutational structure, and a gene regulatory network model, mimicking the mechanisms of gene expression regulation. Simulations showed that correlated mutational effects evolved in the three genetic architectures as a response to correlated selection, but the response in gene networks was specific. The intensity of gene co-expression was mostly explained by the regulatory distance between genes (largest correlations being associated to genes directly interacting with each other), and the sign of co-expression was associated with the nature of the regulation (transcription activation or inhibition). These results concur to the idea that gene network topologies could partly reflect past selection patterns on gene expression.

Kernelized multiview signed graph learning for single-cell RNA sequencing data

BackgroundCharacterizing the topology of gene regulatory networks (GRNs) is a fundamental problem in systems biology. The advent of single cell technologies has made it possible to construct GRNs at finer resolutions than bulk and microarray datasets. However, cellular heterogeneity and sparsity of the single cell datasets render void the application of regular Gaussian assumptions for constructing GRNs. Additionally, most GRN reconstruction approaches estimate a single network for the entire data. This could cause potential loss of information when single cell datasets are generated from multiple treatment conditions/disease states.ResultsTo better characterize single cell GRNs under different but related conditions, we propose the joint estimation of multiple networks using multiple signed graph learning (scMSGL). The proposed method is based on recently developed graph signal processing (GSP) based graph learning, where GRNs and gene expressions are modeled as signed graphs and graph signals, respectively. scMSGL learns multiple GRNs by optimizing the total variation of gene expressions with respect to GRNs while ensuring that the learned GRNs are similar to each other through regularization with respect to a learned signed consensus graph. We further kernelize scMSGL with the kernel selected to suit the structure of single cell data.ConclusionsscMSGL is shown to have superior performance over existing state of the art methods in GRN recovery on simulated datasets. Furthermore, scMSGL successfully identifies well-established regulators in a mouse embryonic stem cell differentiation study and a cancer clinical study of medulloblastoma.

A data-driven optimization method for coarse-graining gene regulatory networks

One major challenge in systems biology is to understand how various genes in a gene regulatory network (GRN) collectively perform their functions and control network dynamics. This task becomes extremely hard to tackle in the case of large networks with hundreds of genes and edges, many of which have redundant regulatory roles and functions. The existing methods for model reduction usually require the detailed mathematical description of dynamical systems and their corresponding kinetic parameters, which are often not available. Here, we present a data-driven method for coarse-graining large GRNs, named SacoGraci, using ensemble-based mathematical modeling, dimensionality reduction, and gene circuit optimization by Markov Chain Monte Carlo methods. SacoGraci requires network topology as the only input and is robust against errors in GRNs. We benchmark and demonstrate its usage with synthetic, literature-based, and bioinformatics-derived GRNs. We hope SacoGraci will enhance our ability to model the gene regulation of complex biological systems.

ScSGL: kernelized signed graph learning for single-cell gene regulatory network inference.

Elucidating the topology of gene regulatory networks (GRNs) from large single-cell RNA sequencing datasets, while effectively capturing its inherent cell-cycle heterogeneity and dropouts, is currently one of the most pressing problems in computational systems biology. Recently, graph learning (GL) approaches based on graph signal processing have been developed to infer graph topology from signals defined on graphs. However, existing GL methods are not suitable for learning signed graphs, a characteristic feature of GRNs, which are capable of accounting for both activating and inhibitory relationships in the gene network. They are also incapable of handling high proportion of zero values present in the single cell datasets. To this end, we propose a novel signed GL approach, scSGL, that learns GRNs based on the assumption of smoothness and non-smoothness of gene expressions over activating and inhibitory edges, respectively. scSGL is then extended with kernels to account for non-linearity of co-expression and for effective handling of highly occurring zero values. The proposed approach is formulated as a non-convex optimization problem and solved using an efficient ADMM framework. Performance assessment using simulated datasets demonstrates the superior performance of kernelized scSGL over existing state of the art methods in GRN recovery. The performance of scSGL is further investigated using human and mouse embryonic datasets. The scSGL code and analysis scripts are available on https://github.com/Single-Cell-Graph-Learning/scSGL. Supplementary data are available at Bioinformatics online.

Reconstruction of dynamic regulatory networks reveals signaling-induced topology changes associated with germ layer specification.

SummaryElucidating regulatory relationships between transcription factors (TFs) and target genes is fundamental to understanding how cells control their identity and behavior. Unfortunately, existing computational gene regulatory network (GRN) reconstruction methods are imprecise, computationally burdensome, and fail to reveal dynamic regulatory topologies. Here, we present Epoch, a reconstruction tool that uses single-cell transcriptomics to accurately infer dynamic networks. We apply Epoch to identify the dynamic networks underpinning directed differentiation of mouse embryonic stem cells (ESCs) guided by multiple signaling pathways, and we demonstrate that modulating these pathways drives topological changes that bias cell fate potential. We also find that Peg3 rewires the pluripotency network to favor mesoderm specification. By integrating signaling pathways with GRNs, we trace how Wnt activation and PI3K suppression govern mesoderm and endoderm specification, respectively. Finally, we identify regulatory circuits of patterning and axis formation that distinguish in vitro and in vivo mesoderm specification.

How to Use the TDCor Algorithm to Infer Gene Regulatory Networks from Time Series Transcriptomic Data.

Over the last few decades, many genes have been functionally characterized and shown to be involved in various metabolic, developmental, and signaling pathways. However it still remains unclear how all these genes and pathways integrate into a unique regulatory network to coordinate the development and the growth, or the response to the environment. This is why unraveling the topology of gene regulatory networks (GRN) has become central to our understanding of all these processes. The recent advancement of high-throughput methods has provided enormous amount of -omics data. These data can now be exploited for rapid network reconstruction with statistical inference methods. We recently published a new GRN inference algorithm called TDCor which reconstructs GRN from time-series transcriptomic data. The algorithm has been released in the form of an R package. Here, I describe into details how to install and use the package.

An Ensemble Method to Reconstruct Gene Regulatory Networks Based on Multivariate Adaptive Regression Splines.

Gene regulatory networks (GRNs) play a key role in biological processes. However, GRNs are diverse under different biological conditions. Reconstructing gene regulatory networks (GRNs) from gene expression has become an important opportunity and challenge in the past decades. Although there are a lot of existing methods to infer the topology of GRNs, such as mutual information, random forest, and partial least squares, the accuracy is still low due to the noise and high dimension of the expression data. In this paper, we introduce an ensemble Multivariate Adaptive Regression Splines (MARS) based method to reconstruct the directed GRNs from multifactorial gene expression data, called PBMarsNet. PBMarsNet incorporates part mutual information (PMI) to pre-weight the candidate regulatory genes and then uses MARS to detect the nonlinear regulatory links. Moreover, we apply bootstrap to run the MARS multiple times and average the outputs of each MARS as the final score of regulatory links. The results on DREAM4 challenge and DREAM5 challenge datasets show PBMarsNet has a superior performance and generalization over other state-of-the-art methods.

Dynamic Modeling of Transcriptional Gene Regulatory Networks.

Diverse cellular phenotypes are determined by groups of transcription factors (TFs) and other regulators that influence each others' gene expression, forming transcriptional gene regulatory networks (GRNs). In many biological contexts, especially in development and associated diseases, the expression of the genes in GRNs is not static but evolves in time. Modeling the dynamics of GRN state is an important approach for understanding diverse cellular phenomena such as cell-fate specification, pluripotency and cell-fate reprogramming, oncogenesis, and tissue regeneration. In this protocol, we describe how to model GRNs using a data-driven dynamic modeling methodology, gene circuits. Gene circuits do not require knowledge of the GRN topology and connectivity but instead learn them from training data, making them very general and applicable to diverse biological contexts. We utilize the MATLAB-based gene circuit modeling software Fast Inference of Gene Regulation (FIGR) for training the model on quantitative gene expression data and simulating the GRN. We describe all the steps in the modeling life cycle, from formulating the model, training the model using FIGR, simulating the GRN, to analyzing and interpreting the model output. This protocol highlights these steps with the example of a dynamical model of the gap gene GRN involved in Drosophila segmentation and includes example MATLAB statements for each step.

DNF: A differential network flow method to identify rewiring drivers for gene regulatory networks

Differential network analysis has become an important approach in identifying driver genes in development and disease. However, most studies capture only local features of the underlying gene-regulatory network topology. These approaches are vulnerable to noise and other changes which mask driver-gene activity. Therefore, methods are urgently needed which can separate the impact of true regulatory elements from stochastic changes and downstream effects. We propose the differential network flow (DNF) method to identify key regulators of progression in development or disease. Given the network representation of consecutive biological states, DNF quantifies the essentiality of each node by differences in the distribution of network flow, which are capable of capturing comprehensive topological differences from local to global feature domains. DNF achieves more accurate driver-gene identification than other state-of-the-art methods when applied to four human datasets from The Cancer Genome Atlas and three single-cell RNA-seq datasets of murine neural and hematopoietic differentiation. Furthermore, we predict key regulators of crosstalk between separate networks underlying both neuronal differentiation and the progression of neurodegenerative disease, among which App is predicted as a driver gene of neural stem cell differentiation. Our method is a new approach for quantifying the essentiality of genes across networks of different biological states.

Characterization of DNA Methylation Associated Gene Regulatory Networks During Stomach Cancer Progression.

DNA methylation plays a critical role in tumorigenesis through regulating oncogene activation and tumor suppressor gene silencing. Although extensively analyzed, the implication of DNA methylation in gene regulatory network is less characterized. To address this issue, in this study we performed an integrative analysis on the alteration of DNA methylation patterns and the dynamics of gene regulatory network topology across distinct stages of stomach cancer. We found the global DNA methylation patterns in different stages are generally conserved, whereas some significantly differentially methylated genes were exclusively observed in the early stage of stomach cancer. Integrative analysis of DNA methylation and network topology alteration yielded several genes which have been reported to be involved in the progression of stomach cancer, such as IGF2, ERBB2, GSTP1, MYH11, TMEM59, and SST. Finally, we demonstrated that inhibition of SST promotes cell proliferation, suggesting that DNA methylation-associated SST suppression possibly contributes to the gastric cancer progression. Taken together, our study suggests the DNA methylation-associated regulatory network analysis could be used for identifying cancer-related genes. This strategy can facilitate the understanding of gene regulatory network in cancer biology and provide a new insight into the study of DNA methylation at system level.

The evolution of gene regulatory networks controlling Arabidopsis thaliana L. trichome development

BackgroundThe variation in structure and function of gene regulatory networks (GRNs) participating in organisms development is a key for understanding species-specific evolutionary strategies. Even the tiniest modification of developmental GRN might result in a substantial change of a complex morphogenetic pattern. Great variety of trichomes and their accessibility makes them a useful model for studying the molecular processes of cell fate determination, cell cycle control and cellular morphogenesis. Nowadays, a large number of genes regulating the morphogenesis of A. thaliana trichomes are described. Here we aimed at a study the evolution of the GRN defining the trichome formation, and evaluation its importance in other developmental processes.ResultsIn study of the evolution of trichomes formation GRN we combined classical phylogenetic analysis with information on the GRN topology and composition in major plants taxa. This approach allowed us to estimate both times of evolutionary emergence of the GRN components which are mainly proteins, and the relative rate of their molecular evolution. Various simplifications of protein structure (based on the position of amino acid residues in protein globula, secondary structure type, and structural disorder) allowed us to demonstrate the evolutionary associations between changes in protein globules and speciations/duplications events. We discussed their potential involvement in protein-protein interactions and GRN function.ConclusionsWe hypothesize that the divergence and/or the specialization of the trichome-forming GRN is linked to the emergence of plant taxa. Information about the structural targets of the protein evolution in the GRN may predict switching points in gene networks functioning in course of evolution. We also propose a list of candidate genes responsible for the development of trichomes in a wide range of plant species.

Tree-Based Learning of Regulatory Network Topologies and Dynamics with Jump3.

Inference of gene regulatory networks (GRNs) from time series data is a well-established field in computational systems biology. Most approaches can be broadly divided in two families: model-based and model-free methods. These two families are highly complementary: model-based methods seek to identify a formal mathematical model of the system. They thus have transparent and interpretable semantics but rely on strong assumptions and are rather computationally intensive. On the other hand, model-free methods have typically good scalability. Since they are not based on any parametric model, they are more flexible than model-based methods, but also less interpretable.In this chapter, we describe Jump3, a hybrid approach that bridges the gap between model-free and model-based methods. Jump3 uses a formal stochastic differential equation to model each gene expression but reconstructs the GRN topology with a nonparametric method based on decision trees. We briefly review the theoretical and algorithmic foundations of Jump3, and then proceed to provide a step-by-step tutorial of the associated software usage.

Transcriptional Profiling of Laser Captured Neurons in the Dorsal Motor Nucleus of the vagus in Response to Ischemic Heart Failure

The dorsal motor nucleus of the vagus (DMV) has been shown to exert a protective effect against damage from myocardial infarction (MI) that may be crucial in preventing subsequent heart failure. A time series of transcriptional changes in 1000 laser captured neurons were studied in a rat model of MI. The results show a significant response after 1 week involving Npff and GABA receptor gene expression changes in neurons. After 3 weeks, there were changes in somatostatin gene expression not seen at 1 week, suggesting increased inhibitory interneuron activity. The was an overall constraining of the diversity of transnational phenotype in response to ischemia both at one week and three weeks. Several of these responses showed left/right asymmetry, with a bias toward changes on the right. Gene regulatory network topology was also significantly altered, showing the potential for the orexin-2 receptor to be a novel regulatory hub. This represents the first time single neurons in the DMV have been assayed in response to ischemic heart failure and the results suggest new ways in which the DMV modulates protective effects at the heart. Support or Funding Information Funding Sources: T32 AA007463-28; U01 HL133360; RO1 HL111621-01A1; OT2 - SPARC The DMV mediates a cardioprotective effect via multiple outputs not only through direct connections with the heart, but also via connections to secretory endocrine cells in gut. This abstract is from the Experimental Biology 2018 Meeting. There is no full text article associated with this abstract published in The FASEB Journal.

Understanding Robust Adaptation Dynamics of Gene Regulatory Network.

Robust adaptation is a critical attribute for gene regulatory network (GRN), understanding the relationship between adaptation and the GRN topology, and corresponding parameters is a challenging issue. The work in this paper includes: first, seven constraint multiobjective optimization algorithms are used to find sufficient solutions to get more reliable statistic rules. Meanwhile, the algorithms are compared to facilitate the future algorithm selection; second, a fuzzy c-mean algorithm is used to analyze solutions and to classify the solutions into different groups; third, the histogram analysis for all satisfactory solutions shows the preferred parameter range, i.e., parameter motif. The contributions of this paper includes: 1) Two new adaptation indices i.e., peak time and settle down time, are proposed for the first time to give more accurate description of the robust adaptation. Our conclusion is that some solutions even with satisfactory sensitivity and precision are not practically of robust adaptation because of too long time needed. 2) The relationship between topology, parameter set, and robust adaptation of GRN is discovered in the sense of both preferred topology and parameter motif. Our conclusion is that the robust adaptation depends more on the GRN topology than the model parameter set in two feasible topologies.

Genome-wide use of high- and low-affinity Tbrain transcription factor binding sites during echinoderm development

Sea stars and sea urchins are model systems for interrogating the types of deep evolutionary changes that have restructured developmental gene regulatory networks (GRNs). Although cis-regulatory DNA evolution is likely the predominant mechanism of change, it was recently shown that Tbrain, a Tbox transcription factor protein, has evolved a changed preference for a low-affinity, secondary binding motif. The primary, high-affinity motif is conserved. To date, however, no genome-wide comparisons have been performed to provide an unbiased assessment of the evolution of GRNs between these taxa, and no study has attempted to determine the interplay between transcription factor binding motif evolution and GRN topology. The study here measures genome-wide binding of Tbrain orthologs by using ChIP-sequencing and associates these orthologs with putative target genes to assess global function. Targets of both factors are enriched for other regulatory genes, although nonoverlapping sets of functional enrichments in the two datasets suggest a much diverged function. The number of low-affinity binding motifs is significantly depressed in sea urchins compared with sea star, but both motif types are associated with genes from a range of functional categories. Only a small fraction (∼10%) of genes are predicted to be orthologous targets. Collectively, these data indicate that Tbr has evolved significantly different developmental roles in these echinoderms and that the targets and the binding motifs in associated cis-regulatory sequences are dispersed throughout the hierarchy of the GRN, rather than being biased toward terminal process or discrete functional blocks, which suggests extensive evolutionary tinkering.

Gene regulatory network inference using PLS-based methods.

BackgroundInferring the topology of gene regulatory networks (GRNs) from microarray gene expression data has many potential applications, such as identifying candidate drug targets and providing valuable insights into the biological processes. It remains a challenge due to the fact that the data is noisy and high dimensional, and there exists a large number of potential interactions.ResultsWe introduce an ensemble gene regulatory network inference method PLSNET, which decomposes the GRN inference problem with p genes into p subproblems and solves each of the subproblems by using Partial least squares (PLS) based feature selection algorithm. Then, a statistical technique is used to refine the predictions in our method. The proposed method was evaluated on the DREAM4 and DREAM5 benchmark datasets and achieved higher accuracy than the winners of those competitions and other state-of-the-art GRN inference methods.ConclusionsSuperior accuracy achieved on different benchmark datasets, including both in silico and in vivo networks, shows that PLSNET reaches state-of-the-art performance.Electronic supplementary materialThe online version of this article (doi:10.1186/s12859-016-1398-6) contains supplementary material, which is available to authorized users.

Netter: re-ranking gene network inference predictions using structural network properties.

BackgroundMany algorithms have been developed to infer the topology of gene regulatory networks from gene expression data. These methods typically produce a ranking of links between genes with associated confidence scores, after which a certain threshold is chosen to produce the inferred topology. However, the structural properties of the predicted network do not resemble those typical for a gene regulatory network, as most algorithms only take into account connections found in the data and do not include known graph properties in their inference process. This lowers the prediction accuracy of these methods, limiting their usability in practice.ResultsWe propose a post-processing algorithm which is applicable to any confidence ranking of regulatory interactions obtained from a network inference method which can use, inter alia, graphlets and several graph-invariant properties to re-rank the links into a more accurate prediction. To demonstrate the potential of our approach, we re-rank predictions of six different state-of-the-art algorithms using three simple network properties as optimization criteria and show that Netter can improve the predictions made on both artificially generated data as well as the DREAM4 and DREAM5 benchmarks. Additionally, the DREAM5 E.coli. community prediction inferred from real expression data is further improved. Furthermore, Netter compares favorably to other post-processing algorithms and is not restricted to correlation-like predictions. Lastly, we demonstrate that the performance increase is robust for a wide range of parameter settings. Netter is available at http://bioinformatics.intec.ugent.be.ConclusionsNetwork inference from high-throughput data is a long-standing challenge. In this work, we present Netter, which can further refine network predictions based on a set of user-defined graph properties. Netter is a flexible system which can be applied in unison with any method producing a ranking from omics data. It can be tailored to specific prior knowledge by expert users but can also be applied in general uses cases. Concluding, we believe that Netter is an interesting second step in the network inference process to further increase the quality of prediction.Electronic supplementary materialThe online version of this article (doi:10.1186/s12859-016-0913-0) contains supplementary material, which is available to authorized users.

Prediction of disease-gene-drug relationships following a differential network analysis.

Great efforts are being devoted to get a deeper understanding of disease-related dysregulations, which is central for introducing novel and more effective therapeutics in the clinics. However, most human diseases are highly multifactorial at the molecular level, involving dysregulation of multiple genes and interactions in gene regulatory networks. This issue hinders the elucidation of disease mechanism, including the identification of disease-causing genes and regulatory interactions. Most of current network-based approaches for the study of disease mechanisms do not take into account significant differences in gene regulatory network topology between healthy and disease phenotypes. Moreover, these approaches are not able to efficiently guide database search for connections between drugs, genes and diseases. We propose a differential network-based methodology for identifying candidate target genes and chemical compounds for reverting disease phenotypes. Our method relies on transcriptomics data to reconstruct gene regulatory networks corresponding to healthy and disease states separately. Further, it identifies candidate genes essential for triggering the reversion of the disease phenotype based on network stability determinants underlying differential gene expression. In addition, our method selects and ranks chemical compounds targeting these genes, which could be used as therapeutic interventions for complex diseases.

Abstract B2-40: Many-attractor models of signaling dynamics and heterogeneity in acute myeloid leukemia gene regulatory networks

Abstract Introduction and purpose of the study: Significant progress has been made in biological network reconstruction methods in recent years, with much emphasis placed on investigating the topology of gene regulatory networks (GRNs). While studying a network's topology can provide useful biological information, a dynamical model for how genes and proteins exchange information is needed in order to understand and predict a cell's response to stimuli such as drugs which inhibit the activity of a protein. Two thirds of patients with acute myeloid leukemia (AML) have an unfavorable prognosis. This is in part due to the high tumor cell heterogeneity in AML, which is often the ultimate cause of drug resistance and relapse in patients. Mathematical models of signaling dynamics in AML which account for the heterogeneity of clonally derived cells in a tumor could be valuable new tools for designing effective, original therapies. Heterogeneity can be described in signaling by defining nonlinear models with multiple attractor states. Novel computational methods: We present two mathematical signaling models that encode real gene expression measurements as attractors in a directed AML GRN. Gene expression profiles obtained from RNA-seq in normal progenitor and AML cells data are used to define a set of robust gene expression profiles corresponding to different clonal states. The first model is Boolean, and is based on an asymmetric Hopfield model with multiple memory patterns (Szedlak et al. 2014). The second model uses continuous expression values in which the elements of the GRN interact like oscillators characterized by multi-stability. The equations of motion for the oscillators are defined in such a way that the expression profiles for the RNA-seq data match the multiple steady states of the signaling network. Summary of new data: We examine the sensitivity of normal and cancer attractors to single gene perturbations and to real gene-inhibiting drugs in both models. For the network topology we use a recently developed network that is specific for AML, specifically AML 2.3 (Ong et al. 2014). Attractors are defined using RNA-seq data from AML cells and hematopoietic controls (Macrae et al. 2013). According to the Hopfield model, we found that ELAVL1, GATA1, IRX5, MYOG, RXRA, and TFEB are genes in AML which, when inhibited, strongly destabilize the cancer attractor. In the continuum model, we explicitly included interactions between drugs currently in AML clinical trials and their targets. Focusing on known targets of lenalidomide and sorafenib, we found that FLT1, KDR, and PDGFRB are associated with the strongest sensitivity to perturbations. In addition, we found that besides the direct targets of these drugs (BRAF, RAF1, FLT4, KDR, FLT3, PDGFRB, KIT, FGFR1, RET, FLT1 for sorafenib, and TNF, TNFSF11, CDH5, PTGS2, CRBN for lenalidomide), RPS18, RPS11, RPS3, and TPT1 are also strongly indirectly down-regulated. Conclusions: We developed two novel methods of network signaling that encode gene expression patterns of cells under varying conditions as attractor states. The methods were applied to AML to identify a set of genes whose perturbation is associated with strong deviations from cancer conditions. Citation Format: Anthony D. Szedlak, Giovanni Paternostro, Carlo Piermarocchi. Many-attractor models of signaling dynamics and heterogeneity in acute myeloid leukemia gene regulatory networks. [abstract]. In: Proceedings of the AACR Special Conference on Computational and Systems Biology of Cancer; Feb 8-11 2015; San Francisco, CA. Philadelphia (PA): AACR; Cancer Res 2015;75(22 Suppl 2):Abstract nr B2-40.

Evolving robust gene regulatory networks.

Design and implementation of robust network modules is essential for construction of complex biological systems through hierarchical assembly of ‘parts’ and ‘devices’. The robustness of gene regulatory networks (GRNs) is ascribed chiefly to the underlying topology. The automatic designing capability of GRN topology that can exhibit robust behavior can dramatically change the current practice in synthetic biology. A recent study shows that Darwinian evolution can gradually develop higher topological robustness. Subsequently, this work presents an evolutionary algorithm that simulates natural evolution in silico, for identifying network topologies that are robust to perturbations. We present a Monte Carlo based method for quantifying topological robustness and designed a fitness approximation approach for efficient calculation of topological robustness which is computationally very intensive. The proposed framework was verified using two classic GRN behaviors: oscillation and bistability, although the framework is generalized for evolving other types of responses. The algorithm identified robust GRN architectures which were verified using different analysis and comparison. Analysis of the results also shed light on the relationship among robustness, cooperativity and complexity. This study also shows that nature has already evolved very robust architectures for its crucial systems; hence simulation of this natural process can be very valuable for designing robust biological systems.